function [C,curvelen,xx_int,xx_ext] = get_geometry_ellipsoid(ntot,xxc,a,b,angle)
% ellipsoid
% ntot: number of nodes
% xxc: center of ellipsoid
% a: radius of ellipsoid in x
% b: radius of ellipsoid in y
% angle: angle of rotation in radians
rotMat=[cos(angle) -sin(angle);sin(angle) cos(angle)];
tt       = linspace(0,2*pi*(1 - 1/ntot),ntot);
C        = zeros(6,ntot );
C(1,:)   =  a*cos(tt);
C(2,:)   = -a*sin(tt);
C(3,:)   = -a*cos(tt);
C(4,:)   =  b*sin(tt);
C(5,:)   =  b*cos(tt);
C(6,:)   = -b*sin(tt);

C([1 4],:)=rotMat*C([1 4],:);
C([2 5],:)=rotMat*C([2 5],:);
C([3 6],:)=rotMat*C([3 6],:);

curvelen = pi*(a+b)*(1+(3*((a-b)/(a+b))^2)/(10+sqrt(4-3*((a-b)/(a+b))^2)));

rmin = sqrt(min(C(1,:).^2 + C(4,:).^2));
rmax = sqrt(max(C(1,:).^2 + C(4,:).^2));

%%% Construct interior points.
ttint  = 2*pi*sort(rand(1,3));
xx_int = 0.5*rmin*[cos(ttint); sin(ttint)];

%%% Construct exterior points.
ttext  = 2*pi*sort(rand(1,5));
xx_ext = 1.5*rmax*[cos(ttext); sin(ttext)];

%%% Shift to center at xxc.
C([1,4],:) = C([1,4],:) + xxc*ones(1,size(C,2));
xx_int     = xx_int     + xxc*ones(1,size(xx_int,2));
xx_ext     = xx_ext     + xxc*ones(1,size(xx_ext,2));

return